
###############################################################################
### Did 3G Make Afghan Insurgents Fight More Effectively? 
### Replication files for analyses
### Mehmet Erdem Arslan
###
### Functions - Spatial Regressions  
###############################################################################

require(dplyr)
require(tidyr)

#### Function for linear combination test ####
lct <- function(mod, var1, var2) {
  mat <- mod$vcov
  coef1 <- mod$coefficients[var1]
  coef2 <- mod$coefficients[var2]
  rownames(mat) <- colnames(mat) <- names(mod$coefficients)
  t <- (coef1 - coef2) / sqrt(mat[var2, var2] + mat[var1, var1] - 2*mat[var2, var1])
  return(round(t[[1]], 3))
}

#### Function for printing coefficient tables for splm & splm_ML class objects ####
# (specific for the article, not for general use)

coef.table <- function(models = list(), bin = FALSE) {
  for (i in (1:length(models))) {
    mod <- models[[i]]
    if ("splm" %in% class(mod)) {
      tab <- data.frame(round(summary(mod)$CoefTable, 3))
      if (!("lambda" %in% rownames(tab))) {
        tab <- rbind(tab, data.frame(round(summary(mod)$ARCoefTable, 3)), 
                     data.frame(round(summary(mod)$ErrCompTable, 3)))
      }
    } else {
      tab <- data.frame(round(summary(mod)$coefficients, 3))
    }
    colnames(tab) <- c("est", "se", "t", "p")
    tab <- cbind(vars = rownames(tab), tab)
    tab <- pivot_longer(tab, cols = c(2,3), names_to = "name", values_to = "value")
    tab$star <- case_when(tab$p <= 0.001 ~ "^{***}", 
                          tab$p <= 0.01 & tab$p > 0.001 ~ "^{**}",
                          tab$p <= 0.05 & tab$p > 0.01 ~ "^{*}", 
                          tab$p > 0.05 ~ NA)
    tab[tab$name == "se", ]$star <- NA
    tab$value <- as.character(tab$value)
    tab$value <- ifelse(is.na(tab$star), tab$value, 
                        paste0(tab$value, tab$star))
    tab[tab$name == "se", ]$value <- paste0("(", tab[tab$name == "se", ]$value, ")")
    tab$star <- tab$t <- tab$p <- NULL
    colnames(tab) <- c("vars", "name", paste0("Model ", as.character(i)))
    if (i == 1) {out <- tab}
    if (i != 1) {out <- merge(out, tab, by = c("vars", "name"), all = TRUE)}
  }
  out$vars <- factor(out$vars, levels = c("(Intercept)", "lag.two", "lag.two.bin", 
                                          "lag.three", "lag.three.bin", "lag.troops", 
                                          "log(plm::lag(all, 1) + 1)",
                                          "log(plm::lag(crd, 1) + 1)", 
                                          "log(plm::lag(ied, 1) + 1)", 
                                          "log(plm::lag(alt, 1) + 1)",
                                          "lag.two:lag.troops", "lag.two.bin:lag.troops", 
                                          "lag.troops:lag.three", "lag.troops:lag.three.bin",
                                          "lambda", "Wtwo", "Wtwo.bin", 
                                          "Wthree", "Wthree.bin", "rho"))
  out <- out[order(out$vars), ]
  out$vars <- as.character(out$vars)
  out[out$name == "se", ]$vars <- ""
  out$vars <- case_when(out$vars == "(Intercept)" ~ "Intercept",
                        out$vars == "lag.two" ~ "2G Coverage$_{t-1}$", 
                        out$vars == "lag.two.bin" ~ "2G Coverage$_{t-1}$ (binary)", 
                        out$vars == "lag.three" ~ "3G Coverage$_{t-1}$", 
                        out$vars == "lag.three.bin" ~ "3G Coverage$_{t-1}$ (binary)", 
                        out$vars == "lag.troops" ~ "Troop presence$_{t-1}$", 
                        out$vars == "log(plm::lag(all, 1) + 1)" ~ "\\emph{log} Attacks$_{t-1}$",
                        out$vars == "log(plm::lag(crd, 1) + 1)" ~ "\\emph{log} Attacks$_{t-1}$", 
                        out$vars == "log(plm::lag(ied, 1) + 1)" ~ "\\emph{log} Attacks$_{t-1}$", 
                        out$vars == "log(plm::lag(alt, 1) + 1)" ~ "\\emph{log} Attacks$_{t-1}$", 
                        out$vars == "lag.two:lag.troops" ~ "2G$_{t-1}$ \\times Troops$_{t-1}$", 
                        out$vars == "lag.two.bin:lag.troops" ~ "2G$_{t-1}$ \\times Troops$_{t-1}$", 
                        out$vars == "lag.troops:lag.three" ~ "3G$_{t-1}$ \\times Troops$_{t-1}$", 
                        out$vars == "lag.troops:lag.three.bin" ~ "3G$_{t-1}$ \\times Troops$_{t-1}$",
                        out$vars == "lambda" ~ "$\\rho$ (Spatial lag of DV)", 
                        out$vars == "Wtwo" ~ "$\\theta_{2G}$ (Spatial lag of 2G)", 
                        out$vars == "Wtwo.bin" ~ "$\\theta_{2G}$ (Spatial lag of 2G)", 
                        out$vars == "Wthree" ~ "$\\theta_{3G}$ (Spatial lag of 3G)", 
                        out$vars == "Wthree.bin" ~ "$\\theta_{2G}$ (Spatial lag of 3G)", 
                        out$vars == "rho" ~ "$\\lambda$ (Spatial error parameter)", 
                        out$vars == "" ~ "")
  out$name <- NULL
  rownames(out) <- NULL
  out <- rename(out, Variables = vars)
  out[is.na(out)] <- " "
  if (bin == TRUE) {
    out <- rbind(out, 
                 c("$t_{\\beta_{3G} > \\beta_{2G}}$", 
                   lapply(models[1:3], FUN = lct, var1 = "lag.three.bin", var2 = "lag.two.bin"), 
                   lapply(models[4:6], FUN = lct, var1 = "lag.troops:lag.three.bin", var2 = "lag.two.bin:lag.troops"))) 
  } else {
  out <- rbind(out, 
               c("$t_{\\beta_{3G} > \\beta_{2G}}$", 
                 lapply(models[1:3], FUN = lct, var1 = "lag.three", var2 = "lag.two"), 
                 lapply(models[4:6], FUN = lct, var1 = "lag.troops:lag.three", var2 = "lag.two:lag.troops"))) 
  }
  print(out, justify = "left", right = FALSE, row.names = FALSE)
}


#### Function to calculate Akaike Information Criterion (AIC) for splm class objects ####
aic_splm <- function(model, k = 2) {
  aic <- (-2*model$logLik) + 2*length(model$coefficients)
  out <- list(AIC = aic, logLik = model$logLik)
  return(out)
}
